active_sites=[1 2 4 7 8];
passive_sites=[9];
control_sites=[10];

mean_resp_by_site=nan(83,9);
aaa=[];
for curr_exp = [1:10]
    % % time_point = 1;
    for time_point = 1:7
        cur_act=[];
        for layer = 1:4
            cur_act(end+1:end+size(proj_meta(curr_exp).rd(layer,time_point).act,1),:) = proj_meta(curr_exp).rd(layer,time_point).act;
        end
        % generate binary activity during feedback.
        velM_bi = proj_meta(curr_exp).rd(1,time_point).velM_smoothed;%(1:7500);
        velM_bi = -velM_bi; % make the velocity of forward-running positive.
        velM_bi(velM_bi > 0.001) = 1;
        velM_bi(velM_bi <= 0.001) = 0;
        
        % generate the frames of ps onset.
        ps_ons = find(diff(proj_meta(curr_exp).rd(1,time_point).ps_id)>1);
        % remove early & late ps.
        ps_ons(ps_ons > 7400) = [];
        ps_ons(ps_ons < 6) = [];
        
        rps_nbr = 0; % number of real ps
        rps_id = [];
        for psid = 1:length(ps_ons)
            if sum(velM_bi(ps_ons(psid)-5:ps_ons(psid)+15)) >= 21
                rps_nbr = rps_nbr+1;
                rps_id(rps_nbr) = psid;
            end
        end

        for ind = 1:length(rps_id)
            rps_frame_id = ps_ons(rps_id(ind));
            tmp = mean(cur_act(:,rps_frame_id-7:rps_frame_id+75));
            aaa(curr_exp,time_point,ind,1)=mean(tmp(10:30))-mean(tmp(1:7));
            aaa(curr_exp,time_point,ind,2)=rps_frame_id;
        end
    end
end


aaa(aaa==0)=nan;

figure; hold on
plot(smooth2(squeeze(nanmean(nanmean(aaa(active_sites,[1:4],:,1),1),2)),1))
plot(smooth2(squeeze(nanmean(nanmean(aaa(passive_sites,[1:4],:,1),1),2)),1),'r')



site_id=6;


act_rsp=aaa(active_sites,[1:4],:,1);
act_rsp=reshape(act_rsp,prod(size(act_rsp)),1);
act_fra=aaa(active_sites,[1:4],:,2);
act_fra=reshape(act_fra,prod(size(act_fra)),1);

act_fra=act_fra(~isnan(act_rsp));
act_rsp=act_rsp(~isnan(act_rsp));


pas_rsp=aaa(passive_sites,[1:4],:,1);
pas_rsp=reshape(pas_rsp,prod(size(pas_rsp)),1);
pas_fra=aaa(passive_sites,[1:4],:,2);
pas_fra=reshape(pas_fra,prod(size(pas_fra)),1);

pas_fra=pas_fra(~isnan(pas_rsp));
pas_rsp=pas_rsp(~isnan(pas_rsp));

figure;
hold on
plot(act_fra,act_rsp,'.')
plot(pas_fra,pas_rsp,'r.')
plot([1 7000],polyval(polyfit(act_fra,act_rsp,1),[1 7000]))
plot([1 7000],polyval(polyfit(pas_fra,pas_rsp,1),[1 7000]),'r')




